Μερικές Διαφορικές Εξισώσεις

Εξίσωση θερμότητας ($\partial_{t} u-\sigma \nabla^2 u=0$)

Σε πεπερασμένο χωρίο

Στήσιμο ΜΔΕ


           Clear["Global`*"]
PDE = D[u[x, t], t] - D[u[x, t], {x, 2}] == 0;
u[x, t] = X[x]*T[t]
PDE = D[u[x, t], t] - D[u[x, t], {x, 2}] == 0

Αρχικές συνθήκες


            init = u[x, 0] == f[x]

Συνωριακές συνθήκες


             bound1 = u[0, t] == 0
bound2 = u[L, t] == 0


             PDE[[1]]/(X[x] T[t]) // Apart

Έχουμε


              ODEt = Derivative[1][T][t]/T[t] == λ
ODEx = X''[x]/X[x] == λ

Όπως και στην περίπτωση της εξίσωσης κύματος, οι περιπτώσεις $λ\geq0$ απορρίπτονται. Άρα:


               λ = -k^2
ODEx
ODEt
ODEx = X''[x] + k^2 X[x] == 0
ODEt = T'[t] + k^2 T[t] == 0

Λόγω των συνοριακών συνθηκών έχουμε πάλι:


                DSolve[{ODEx, X[0] == 0, X[L] == 0}, X[x], x]

Έχουμε ιδιοτιμές τις (nπ)/L και ιδιοσυναρτήσεις τις sin((nπ x)/L).


                 k = (n Pi)/L
DSolve[{ODEx, X[0] == 0, X[L] == 0}, X[x], x]
DSolve[ODEt, T[t], t]

Έχουμε ότι η $u_n(x,t)=\exp(-\frac{n^2 \pi^2 t}{L})\sin(\frac{n \pi x}{L})$ ικανοποιεί τη Μ.Δ.Ε. και τις αρχικές συνθήκες. Το ίδιο συμβαίνει και με κάθε γραμμικό συνδυασμό της. Θα ψάξουμε να βρούμε ποιος ικανοποιεί και την αρχική συνθήκη. Θέλουμε: $\sum_{n=1}^{\infty} c_n u_n (x,0)=f (x)\Leftrightarrow \sum_{n=1}^{\infty} c_n \sin(\frac{n \pi x}{L})=f(x)$ Πολλαπλασιάζουμε εσωτερικά τα δύο μέλη της τελευταίας ισότητας και λαμβάνουμε υπ' όψιν την ορθογωνιότητα των ημιτόνων. Έτσι:


                    c[n_] := Assuming[Element[n, Integers], Integrate[f[x]*Sin[(n π x)/L], {x, 0, L}]/Integrate[Sin[(n π x)/L]^2, {x, 0, L}]]
c[n]

Έτσι, η προσέγγιση της λύσης με $n_0$ όρους είναι η:


                     un[x_, t_, n_] := c[n] E^(-((n^2 π^2 t)/L^2)) Sin[(n π x)/L]
uApprox[x_, t_, n0_] := Sum[un[x, t, n], {n, 1, n0}]
uApprox[x, t, 4]

Ας δούμε τι πετύχαμε.


                      L = 2 Pi;
f[x_] := x^2


                      uApprox[x, t, 4]


                      Table[Plot3D[Evaluate[uApprox[x, t, n]], {x, 0, L}, {t, 0, 10}], {n,  2, 6}]


                      Plot3D[Evaluate[uApprox[x, t, 20]], {x, 0, L}, {t, 0, 10}, AxesLabel -> {"x","t"}]